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ABSTRACT 

We investigate the large-scale evolution of a relativistic magnetic reconnection 
in an electron-positron pair plasma by a relativistic two-fluid magnetohydrody- 
namic (MHD) code. We introduce an interspecies friction force as an effective 
resistivity to dissipate magnetic fields. We demonstrate that magnetic reconnec- 
tion successfully occurs in our two-fluid system, and that it involves Petschek-type 
bifurcated current layers in a later stage. We further observe a quasi-steady evo- 
lution thanks to an open boundary condition, and find that the Petschek-type 
structure is stable over the long time period. Simulation results and theoreti- 
cal analyses exhibit that the Petschek outflow channel becomes narrower when 
the reconnection inflow contains more magnetic energy, as previously claimed. 
Meanwhile, we flnd that the reconnection rate goes up to ~1 in extreme cases, 
which is faster than previously thought. The role of the resistivity, implications 
for reconnection models in the magnetically dominated limit, and relevance to 
kinetic reconnection works are discussed. 

Subject headings: magnetic flelds — relativity — Magnetohydrodynamics: MHD 
— plasmas 



1. INTRODUCTION 



Magnetic reconnection in coUisionless or collisional plasmas is the driver of explosive 
events in space and astroplasmas. By breaking the magnetic fleld topology, it rapidly re- 
leases the magnetic energy into plasma kinetic energy in a short timescale, and therefore it 
explains particle acceleration or bursty emission signatures in these sites. On the Sun, it is 
widely recognized that magnetic reconnection drives solar flare or coronal mass ejections (see 
Aschwanden (2006) for review). Theoretical models have long been established ( Sweet|[l958 



Parker 1957 Petschek 1964), and a series of MHD simulations make a significant success to 



understand fiare-type events (e.g., Chen & Shibata (2000); Yokoyama & Shibata (2001)). 
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Magnetic reconnection is considered in a wide variety of high-energy astrophysical con- 
texts too. For example, soft gamma repeaters (SGRs) and anomalous X-ray pulsars (AXPs) 
are now best described by a "magnetar" model (Duncan & Thompson 1992 [Woods & 



Thompson 2006), a neutron star with superstrong magnetic fields up to 10 — 10 G. In 



analogy to the Sun, flares on and around the magnetar (Thompson & Duncan 1995, 2001 



Lyutikov 2003 , 2006 ) are considered as driving mechanism of bursty events, in relativistic 



electron-positron environments. Such flares, or magnetic reconnection events, should be 
strongly influenced by the relativistic effects, because the ultra strong magnetic field boosts 
the Alfven speed up to the light speed. 

The pulsar environments are also influenced by relativistic plasmas and the strong mag- 
netic fields (~ lO^^G) of the neutron star. Recent time-dependent simulations of pulsar 



magnetospheres (Komissarov 2006 Bucciantini et al. 2006 Spitkovsky 2006) suggested that 



the magnetic reconnection near the Y point, where the outmost closed field lines intersect the 
equatorial current sheet, is of critically importance, while these models cannot deal with local 
reconnection physics. Outside the magnetosphere, reconnection processes in the "striped" 
current sheets are considered to dissipate magnetic energy inside the relativistic plasma out- 



flow (pulsar winds; Michel (1982 1994); Coroniti (1990); Lyubarsky fc Kirk (2001); Kirk 



& Skjaeraasen (2003)) and its termination shock (Lyubarsky 2003). Furthermore, active 



galactic nuclei ( |di Matteo] 1998 Birk et al. 2001), extragalactic jets (Lesch & Birk 1998), 



gamma-ray burst (GRB) outflows (Drenkhahn 2002; Drenkhahn & Spruit 2002), and po- 



tentially the black hole ergosphere (Koide & Aral 2008) may be influenced by the magnetic 



reconnection in the relativistic regime. Indeed, there is a high demand for modeling the 
magnetic reconnection in these relativistic environments. 

However, the relativistic theory of a magnetic reconnection is not yet well established. 



Blackman & Field ( 1994 ) extended the steady state reconnection models into the relativistic 



regime, based on a relativistic extension of Ohm's law (Blackman & Field 1993). Assum- 



ing uniform proper density, they argued that the Lorentz boost may enhance the energy 
conversion rate both in Sweet-Parker and in Petschek reconnections. In the Sweet-Parker 



regime, Lyutikov & Uzdensky (2003) further examined this idea and claimed that reconnec- 



tion outflow may be super- Alfvenic. On the other hand, Lyubarsky (2005) pointed out that 



the reconnection will not be fast because the relativistic gas pressure increases the outflow 
inertia. Recently, the authors discussed a two-fluid description and we showed that the in- 



compressibility assumption is invalid for relativistic outflow (Zenitani & Hesse 2008b). In 



the Petschek regime, in which the reconnection involves a bifurcated slow-shock structure. 



Lyubarsky (2005) argued that the reconnection would not be an efficient energy converter 



because the slow-shock angle becomes narrower. 
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Meanwhile, there has been a remarkable progress on the kinetic-scale behaviors of rela- 



tivistic magnetic reconnection, by self-consistent particle-in-cell (PIC) simulations. Zenitani 



& Hoshino (2001) demonstrated that powerful DC acceleration occurs around the reconnect- 



ing X-type region. This and the relevant particle acceleration generate nonthermal plasma 



distributions on a larger scale (Jaroschek et al. 2004; Zenitani & Hoshino 2007 Bessho & 



Bhattacharjee 2007 Karlicky 2008), and particle acceleration may be enhanced in a com- 
pressed pulsar-wind configuration ( Lyubarsky fc Liverts||2008 ). In the orthogonal plane, the 
current-driven drift kink instabilities are of importance (Zenitani & Hoshino 2005a 2007), 



because they grow faster and may interfere with the magnetic reconnection. Due to a wide 
variety of such plasma instabilities the reconnection current sheet exhibits complex evolution 



in three dimensions (Jaroschek et al. 2004; Zenitani & Hoshino 2005b, 2008). Furthermore 



it was recently pointed out that kinetic effects are important not only in the critical re- 



connecting region (Hesse & Zenitani 2007), but also in the reconnection outflow region as 



an anisotropy-driven Weibel-type instability (Zenitani & Hesse 2008a). However, these PIC 



simulations typically deal with the spatial domain of several hundreds of the plasma inertial 
length {c/ujp) in two or three dimensions. The large-scale evolution of relativistic reconnec- 
tion systems is still an open problem. 

In order to study large-scale properties of a relativistic magnetic reconnection beyond 
these kinetic scales, and in order to investigate larger scale astrophysical problems which 
contain relativistic magnetic reconnection such as magnetar flares and global pulsar mag- 
netospheres, we need a relativistic extension of magnetohydrodynamic (MHD) codes (see 



Marti & Miiller (2003) for review). However, relativistic hydrodynamic codes are difficult 
to develop, because of the complexity of the equation system. In particular, these codes 
typically use an inverse transformation from the conserved variables in the lab frame to the 
primitive variables in the proper frame. This can be calculated by solving quartic equations, 
or by using iterative methods (e.g. Duncan & Hughes (1994)). Such inverse conversion is 



further complicated in the ideal MHD cases dKoide et al.|[l996| |Komissarov|[l999| [PeT Zianna 



et al. 


2003 


Noble et al. 


2006 



progress both in relativistic magnetohydrodynamic (RMHD) codes and in general relativistic 



magnetohydrodynamic (GRMHD) codes (Koide et al. 1999 Gammie et al. 2003 Mizuno et 



al. 2006). 



To deal with the magnetic reconnection problems, one has to incorporate "resistive" 
effects into the RMHD equations. Otherwise, only the numerical resistivity plays a role to 
dissipate magnetic fields. The first resistive RMHD work was done by Watanabe & Yokoyama 



(2006), by using a spatially limited resistivity. Although their system size is very small (416 



X 200), they successfully presented a Petschek-type reconnection in a mildly relativistic 



regime. Komissarov (2007) also developed the upwind scheme for resistive RMHD, which 
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may be applicable to the reconnection problem. These resistive RMHD studies are based on 



a simple form of time-stationary Ohm's law (Blackman & Field 1993) 



In the present paper, we investigate large-scale properties of a relativistic magnetic re- 
connection in an electron-positron pair plasma by means of two-fluid RMHD simulations. In 
contrast to the conventional RMHD models, we introduce a relativistic two-fluid approxima- 
tion for the first time to our knowledge, so that we can describe the physics in more detail. 
An interspecies friction term is introduced in the momentum equations, which works as an 
effective resistivity. By using a spatially limited resistivity profile, we successfully reproduce 
a magnetic reconnection. We also note that we carry out larger scale simulations, directly 
solving equations to restore the primitive variables. 

This paper is organized as follows. In Section 2, we describe our simulation model. 
Mathematical procedures are also presented in the appendix chapters. In Section 3, we 
overview the system evolution in detail, and present parameter dependences. Especially, we 
analyze the structure of bifurcated Petschek-type current layers in depth. We also demon- 
strate that the system evolution highly depends on the resistivity model. In Section 4, we 
discuss the characteristics of the two-fluid approach and implications for the reconnection in 
the magnetically dominated limit. The last section Section 5 contains the summary. 



2. SIMULATION MODEL 

We employ a relativistic two-fluid model of electrons and positrons. The electron motion 
and positron motion are considered separately. The continuity equation, the momentum 
equation, and the energy equation of relativistic positron fluid, and Maxwell equations are as 
follows. In addition, we introduced an interspecies friction term to the momentum equation, 
which is proportional to the relative motion of electrons and positrons. 

= -Qfpnp = -V-[npUp) (1) 

+lpnpqp{E + — X B)- TfrNpNe{Vp - Ve) (2) 



c 



dKp _ d 



(l^Wp -pp- NpUic^^ 
7 ■ ('jpWpUp — npinc^Up 
-cV X E (4) 



dt dt\ 

= -V ■ {jpWpUp - npinc^Up) + -fpUpqp^Vp ■ E) (3) 

dB 

~dt 
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— = cV X B - 47r ^ qsUsUs (5) 

s=p,e 

In these equations, the subscript s denotes the species ("p" for positrons, and "e" for elec- 
trons), A'^ is the lab-frame density, 7 is the Lorentz factor, n is the proper density, u — 
is the fluid 4-velocity, m is the momentum density, w is the specific enthalpy, 5ij is the 
Kronecker delta, p is the proper isotropic pressure, — —q^ is the positron/electron charge, 
Tfr is the coefficient for an inter-species friction, and K is the kinetic energy density (energy 
density without the rest mass energy). The enthalpy w is defined in the following way: 

w = e -\- p = nmc? -\- \r / (V — l)]p = hnmc? ^ (6) 

where e is the internal energy, F = 4/3 is the specific heat, and h is the dimensionless specific 
enthalpy. 

We solve the equations by using modified Lax-Wendroff scheme. To restore the primitive 
variables 7, u) from the conservative variables m, and K, we use the following 
quartic relation for u — \u\/c 



f(u) = G'{S' - M')u* -2GMDu 



+ 



G'£' - D' -2GM\G -I) 



-2DM{G - l)u - (G - ifM^ = (7) 

where D — Nmc^, M — \m\c, S — K + D, and G — r/(r — 1). We algebraically solve this 
equation by decomposing the quartic equation into the product of two quadratic equations. 
See Appendices A and B for details. We stop the simulation when we find multiple possible 
solutions or when the solution is physically invalid (e.g., negative density). We added small 
artificial viscosity to the code, which works when the fluid 4-velocity has a strong shear so 
that it reduces a numerical oscillation near discontinuities. 

We study the system evolution in the two-dimensional x-z plane. We choose the follow- 
ing relativistic Harris model as an initial configuration: 

B = Botanh{z/L)x (8) 

E = Veffj (9) 

j = 2qpnoUo cosh~^ {z / L) y = ^ qsnsUs{z) (10) 

s=p,e 

Ug = no cosh.~'^ (z / L) + (11) 
Ps = po cosh'^ {z / L) + Pin (12) 
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where L is the typical half-thickness of the current sheet. In the electric field, r/e// = {.Tfr/lp) 
is an effective resistivity, and Mq stands for the initial positron drift to carry the current. 
We also consider uniform background plasmas whose density and pressure are and pin, 
respectively. In this work, the plasma pressure in the Harris sheet is set to po = nomc^. The 
background pressure is set to pin = riinmc^ unless stated otherwise. 

In the case of two-dimensional antiparallel reconnection, we already know that positron 
motion and electron motion are the same in the x-z plane and the opposite in the y-direction. 
Therefore, we assume the following symmetric motion Upx = Uex, Upy = —Uey, Upz = Uez, 
Hp = He, Pp = Pe SO that wc cau reduce the computational cost. Consequently, the current 
has only the ^/-component and = Jz = 0, and we can neglect three components of the 
electromagnetic field, = = By = 0. The assumption also justifies that we do not 
consider the interspecies energy transfer in equation |3| because we assume such a symmetric 
model. In PIC simulations, one characteristic process to generate the charge separation is the 



Weibel instability (Zenitani & Hesse 2008a), driven by an anisotropy in plasma distribution 
function; however, such a small-scale kinetic effect is out of scope of this fluid paper. In a 
MHD-scale, charge neutrality is plausible. By assumption, we do not need to deal with the 
Poisson equations V ■ E = 47r ^^(gs7s?T,s) in this system. 

We introduce a spatially localized resistivity by controlling the interspecies friction force. 
Its profile is set in the following way: 

Tfr = T-o + ri cosh-'[Vx2 + zV(2L)], (13) 

where the background value tq is equivalent to the Reynolds number S = 3000, and localized 
value Ti is equivalent to = 30. In addition, a magnetic field perturbation is added to the 
initial model to quickly trigger a magnetic reconnection. It is defined by the following vector 
potential: 

6 Ay = 2LBi exp[-{x^ + z'^)/{2L)% (14) 

where Bi = 0.03i?o is the typical peak amplitude of the perturbed field. 

The boundaries are located at x = and z = ±Lz, and the reconnection is considered 
around the origin. Boundary conditions for the fluid properties, the electric field, and the 
tangential magnetic field are set to open: d/dx = at the x-boundaries (outflow boundaries) 
and d/dz = at the 2;-boundaries (inflow boundaries). The normal component of the 
magnetic field is set so that it satisfies V ■ -B = at the boundaries. The system size 
{2Lx X 2Lz) is presented in Table [l] in the unit of L. The thickness is typically resolved by 20 
grids {L = 20 Ag). This grid size is selected so that it is comparable to the kinetic scale of a 



typical gyroradius Ag ~ {mc / qpBo) = 0.05L (Equations (15) and (16) in Zenitani & Hoshino 



(2007)). In this equilibrium, the electron inertia length is c/ujp = [mc^/(47r7^no5p)]^^^ 
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Uz[nQmc^ /pQ\^f'^ = 0.1, based on the reference density Uq. The time step is set to At = 
0.2(Ag/c) = 0.01 Tc, where Tc = L/c is the hght transit time. It is sufficiently small, 
{qpBo/hmc)At ~ 0.04 ^ 2n, with respect to the fluid bulk motion. 

Our code is originally developed from the CANS code, a collection of hydrodynamic and 
MHD codes, which has been extensively used in Japanese solar and astrophysical community. 
The code is massively parallelized by MPI. 

We carry out various simulation runs with different parameters. The list of simulation 
runs is presented in Table [Tj The parameter am is the magnetization parameter, which 
stands for the ratio of the magnetic energy flow to the rest mass energy flow, 

' - (15) 



47rm{2'y'^n)c'^ 

Another parameter is the exact ratio of the magnetic energy flow to the plasma energy 
flow, which contains relativistic pressure effect 

The Alfven speed ca in the relativistic regime can be written as follows: 



The subscript in {am,in, o'e,in and CA,in) stands for the upstream values, based on the initial 
inflow properties (e.g. nin,Pin)- Later we often use ae^in as a measure of the upstream energy 
composition. In Table [T| run U3 employs the uniform resistivity model without the ti term 
in equation 13 Runs S3, M3 and XL3 are done in different resolutions. 

Before visiting the simulation results, let us clarify the role of a newly introduced friction 
term. From the positron momentum equation (Equation |2|, we obtain the following relation: 

E + ^xS 

c 

Upnip^Up ■ V)hpUp + niphpUplV ■ (upUp)] 



""^pTlpClp 

+^Pp + -T^lp'mpnphpUp 



TfrNpN, 



Vp - Ve 



nip 



'ypTipQp 

f 9 „\, I „ TfrNpNf, , , 
+ ■ ^ K'^V ^Vv + (^P - *^e)- 

\dt J IpUpQp IpUpqp 
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We consider Ohm's law in the y direction. Dropping d/ dy and considering symmetric electron 
properties, we obtain 

-d{hpUpy) d{hpUpy) 



Ey + 



,Vr. 



X B) 



rup 

+Vp^ 



dt 
d{hpUpy 
dz 



+ v. 



px 



dx 

+ Veffjy 



(18) 



Thus, the fluid inertial effect, the momentum advection, and the interspecies friction term 
work as an effective resistivity. In our two-fluid model, this interspecies resistivity plays an 
essential role to sustain the magnetic reconnection. Around the reconnecting X-point, the 
Lorentz term is negligible because S ~ and Vx,Vz ~ 0, the advection terms usually vanish 
by symmetry, and the inertial terms do not work in the quasi-steady condition {d/dt ~ 0). 
Therefore, the interspecies resistivity sustains the reconnection electric field Ey ~ f]effjy 
Note that the reconnection cannot go on without the reconnection electric field Ey. We do 
not assume any specific mechanism as the interspecies friction term. In a collisional regime, 
it should be equivalent to the collisional term; however, we do not know the true form of the 
relativistic collisional term, which often relies on empirical functions (e.g.. Section 7 in Clare 



& Strottman (1986)). In a coUisionless regime, it is known that the off-diagonal part of the 



pressure tensor sustains the reconnection electric field ( Hesse fc Zenitanl||2007 ) in the kinetic 
simulations. Although its physical meaning is not yet well established, the off-diagonal part 
of the pressure tensor contains several kinetic effects such as the escaping convection of the 
accelerating particles, or the inertial effect of thermal plasma populations. The purpose of 
the interspecies resistivity is to represent these kinetic effects in the fluid approximation, for 
the purpose of larger scale modeling. 



3. SIMULATION RESULTS 

3.1. Evolution overview 

In this section, we overview the system evolution of our reference run (run L3) in detail. 
Due to the trigger field, magnetic reconnection occurs around the center of the simulation 
domain. Plasma outflows start to travel into the ±x directions from the center, while inflows 
come from the ±z directions. The panels in Figure [T] show various physical properties at 
t/Tc = 75 in the normalized unit: the plasma proper density n, the plasma 4- velocity Ur^, 
the electric current jy = 2qpnUy, and the reconnection electric field Ey. Since reconnection 
outflows eject a lot of plasmas, we see dense plasma islands (plasmoids) around x/L ^ ±25- 
30 (Figure [T]a). The reconnection outflow jets become very fast, up to ~ 3.28c (Figure 
[Tfc). The Ux profile shows a characteristic crab claw structure in the plasmoid region, because 
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Fig. 1. — Snapshots of run L3 at I/tc = 75 in the x-z two-dimensional plane, (a) The 
plasma proper density n/no, (6) the x-component of the 4- velocity of the plasma flow Ux/c, 
(c) the out-of-plane electric current jy/jo, and (d) reconnection electric field Ey/Bo. The 
solid lines show magnetic field lines. 
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the dense current sheet plasmas exist around the neutral plane (2; ~ 0). We also see weak 
reverse flows around xjL ^ ±20 after plasmoid passing. In the reconnecting region, there 
is a thin central current layer, and its peak current is 2.5 times larger than the initial state 
(Figure [T]c). Similar enhancement is often seen in classical nonrelativistic models. There 
are rather complicated current structures inside the plasmoids. At this stage the out-of- 
plane electric fleld (or the reconnection electric fleld) is well developed (Figure [iji). The 
typical amplitude is EyjE^ ~ 0.1 over the reconnection region. In addition, the electric 
field is enhanced EyjB^ ~ 0.6 around xjL ^ ±20-25, where reconnected magnetic flux 
is accumulated. The energy and momentum of these enhanced fields are converted to those 
of the downstream plasmas. We also note that this pileup region plays an interesting role 



as particle accelerator (Jaroschek et al. 2004 Zenitani & Hoshino 2007). In general, the 



magnetic topology, electric field properties, and spatial distribution of plasma properties are 
sufficiently consistent with previous reconnection studies by PIC or MHD simulations. The 
system evolution is similar to the Sweet-Parker reconnection which features a single current 
sheet, although the reconnection grows fast. 

After the initial phase, the reconnection continues and plasmoids travel into the ±x- 
directions. The top two panels in Figure [2] show late-time snapshots at t/xc = 200. At this 
stage, plasmoids start to reach the outfiow boundaries, as we see in the profile (Fig. [2]a). 
Note that the entire domain is presented in the x direction. The fastest fiows ~ 3.5c 
are found at xjL ~ ±90 along the outfiow line, where the outfiow channels are connected 
to the plasmoids. An important feature is found in the electric current profile (Figure 
|2^). From the central X-type region to the downstream region, the current layers are now 
bifurcated. The bifurcation starts around x/L = ±40 at I/tc = 100-125. We think these 
current layers are a signature of the Petschek-type steady reconnection, which enables faster 



energy conversion, and we analyze their structure in a later section (see 3.3). Interestingly, 
we see weak "reverse currents" between the two current layers. The current structures inside 
the plasmoids become further complicated, including the interaction with boundaries. 

Since we employ the open boundary condition, plasmoids and reconnection outfiows pass 
through the x-boundaries. Since plasmas and magnetic field lines are continuously supplied 
from the infiow open boundaries at z = ±Lz, the reconnection still continues, and therefore 
the system evolves further. Importantly, the system grows into a steady state reconnection 
structure after the plasmoids have left. The bottom two panels in Figure |2] show the snapshots 
of a very late stage at t/Tc = 400. Now the outfiow channels (Figure |2]c) between Petschek- 
type current layers (Figure |2](i) are found all over the x direction. The distance between the 
two current layers is ~ 2.5-3L at the outfiow boundary {x/L = = 120). Thus, the slope 
angle of the current layer is very small, compared to a typical slow-shock angle of ^ ~ 0.1 in 
nonrelativistic Petschek reconnection. The magnetic field line structure is very smooth over 
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Fig. 2. — Large-scale snapshots of run L3 in the x-z two-dimensional plane, (a) The x- 
component of the plasma 4- velocity u^/c at t/Tc = 200, (b) the out-of-plane current jy/jo at 
t/Tc = 200, (c) the x-component of the plasma 4-velocity Ux/c at t/Tc = 400, and (d) the 
out-of-plane current jy/jo at t/Tc = 400. The black lines show the magnetic field lines. We 



later discuss the properties along the white line {x/L = 100) in panel (d) in Section 3.3 
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the entire simulation domain. We find that these current layers remain stable for a relatively 
long time. 

Next, we investigate the structures of the outflow region in more detail. Figure |3] presents 
the temporal evolution of physical properties along the outflow line (2; = 0) in the normalized 
units. We compare the three nonsteady stages (t/Tc = 75, 150,200) in black lines, and the 
late-time steady stage (t/rc = 400) in red. The vertical magnetic field (Figure [3]a) is a 
reconnected component of magnetic field lines. It is zero at the X-point, and it remains at 
constant level of ~ O.I-Bq inside the outflow channel. Strong peaks are the pileup regions, 
where the reconnected field lines are piled up in front of the dense plasmas. As discussed, 
the electric fields are also enhanced there. Such a powerful magnetic pileup and the relevant 
motional electric fields are signatures of fast magnetic reconnection. The pileup is so strong 
that several discontinuities appear near the pileup regions. For example, at the upstream 
side of the pileup region, the outflow speed Ux becomes very fast but it suddenly drops 
(Figure [3^). On the other hand, there is a strong jump in Bz at the downstream side of the 
pileup region, although the velocity jump is not so clear. We think they are the tangential 
discontinuity or a weak shock (the downstream one) and the relevant reverse fast shock (the 
upstream one). In the later stages (t/rc > 150), the system starts to suffer from numerical 
noises in the downstream side of the plasmoids, as seen in the velocity profile or in the density 
profile (Figure [3]c?). These noises go away as plasmoids pass through the outflow boundaries. 
Importantly, we find that the out-of-plane 4-velocity Uy is not negligible over the relatively 
large region |x/L| < 20 (Figure [3]c). Since Uy is coupled with in-plane components and u^, 
this immediately implies that the conventional one-fluid MHD approximation breaks down 
and that the two-fluid approximation is essential there. At t/xc = 200, Uy becomes negative 
around xjL ^ 80-90. This stands for the negative current between the Petschek-type current 
layers. The bottom panel (Figure [sje) shows the plasma temperature T = p/nmc^ . It is very 
large at the reconnecting X-point, and also (T ~ 2nmc^) inside the outflow channel. The 
typical Lorentz factors in the outflow region are 7 ~ 2.4±0.2 {t/Tc = 200) and 2.1 ± 0.1 
(t/Tc = 400), They are comparable with an Alfvenic value (1 + ae^in)^^'^ = 2.2. 

It is important that the late-time profiles at t/r^ = 400 (indicated by the red lines in 
Figure [sj are quite similar to the earlier profiles. This tells us that the late-time structure 
(Figures [2]c and [2](i) is a very good prediction of the steady state profile. We still see a 
numerical noise around x/L ~ 60-70. This is because this outflow channel is located in the 
downstream side of the shock-type region. As discussed, the outflow channel is located at 
the downstream side of the two current layers. 

Next, we visit the physical properties along the inflow line (x = 0). From Figure |4]a we 
know that the reconnection starts to consume the antiparallel magnetic field B^, but it goes 
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Fig. 3. — Temporal evolution of physical properties along the outflow line {z = 0) in run L3. 
(a) The vertical magnetic fleld Bz/Bq, (b) the outflow component of the positron 4- velocity 
u^/c, (c) the out-of-plane component of the positron 4-velocity Uy/c, (d) the normalized 
plasma number density 'jn/rio, and (e) the normalized plasma temperature p/{nmc^). 
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Fig. 4. — Temporal evolution of physical properties along the inflow line {x — 0) in run L3. 
(a) The antiparallel magnetic field B^/Bq, (b) the reconnection electric field Ey/Bg, and (c) 
the infiow positron velocity Vz/c. 
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down to an asymptotic level of 5^. ~ ±0.74i?o due to the open boundary condition. The 
field reversal is localized in the narrow region around z ~ 0. The reconnection electric field 
Ey grows as the system evolves (Figure [4^). At the later stages, it becomes constant over the 
simulation domain. This tells us that our open-boundary condition works excellently. Also, 
plasma infiow remains at the constant level around the center (Figure |4]c). This us tells that 
the reconnection constantly goes on, consuming outside plasmas and magnetic fields at the 
constant rate. 

The amplitude of the reconnection electric field 



at the X-point is one of the most important parameters in a magnetic reconnection. This 
measures how fast the system transports the magnetic fiux into the X-point, or how fast 
the reconnection consumes the upstream magnetic energy. It is often referred as the "re- 
connection rate" in various normalized form. Following convention, we used the following 
reconnection rate, because reconnection outfiow speed is often approximated by the upstream 
Alfven speed: 



Here the subscript in' denotes the infiow properties measured at z/L = 20. The time 
evolution of r(t) and f{t) is presented in Figure [s] In addition to the reference run L3, two 
other runs M3 and XLS (similar runs with difference resolutions) are overplotted in order to 
check the convergence of the simulation: three are in excellent agreement. The normalized 
rate f{t) is larger than the raw rate r{t), mainly because the inflow magnetic fleld -B^.m' 
decreases over time (Figure |4]a). We see that both the rates remain stable throughout the 
system evolution. Indeed, the normalized rate remains constant: f{t) ~ 0.14. 

The other quantity r* (t) is the time derivative of the accumulated magnetic flux along 
the inflow line 



Because of the discrete sampling time, the calculated value is rather crude, but is useful 
enough to validate the simulation results. In the early stage, both r{t) and r*{t) are in 
excellent agreement. They do not agree after I/tc > 80, because the magnetic flux enters 
from the open inflow boundaries. During t/rc ~ 340-400, r*{t) exhibits strange behavior. 
We conflrmed that this is a boundary effect. Since plasmoid passes through the outflow 
boundaries around t/xc ~ 200-250, perturbation travels from there as a light wave or a 



r{t) = Ey/Bo 



(19) 




(21) 
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fast Alfven wave. The waves from the two outflow boundaries arrived at the center of the 
inflow boundaries. Since two waves carry outward energy flux, the incoming magnetic flux 
temporally slows down, but the system adjusts itself and it goes back to the quasi-steady 
state after I/tc > 400. Note that the final asymptotic value r*{t) = indicates the steady 
evolution. 



CD 

g 
o 

c 
o 
o 
o 

cc 



0.2 



0.15 - 



0.1 - 



0.05 - 



- 



r(t) 




J I I L 



J I I L 



± 



J I I L 



J I I L 







100 200 300 400 
Light Transit Time (L/c) 



500 



600 



Fig. 5. — Time evolution of the reconnection rates in runs L3 {solid lines), M3 {dashed 
lines), and XLS {gray thick lines). The raw reconnection rate r{t) = Ey/ Bq at the X-point, 
the normalized reconnection rate f{t) = cEy/[cA,in'Bx,in'], and the flux consumption rate 
r*{t) are presented. 



3.2. Case studies 

In this section, we compare various simulation runs, focusing on the composition of 
the typical upstream energy flow ae^in- As presented in Table [l| this parameter is mainly 
controlled by the upstream plasma density riin/nQ. The magnetically dominated cases of 
c"e,m 3> 1 ("high-cr" runs) are of strong astrophysical interest, while plasma-dominated cases 
of a^^in < 1 (low-cr runs) can be compared with nonrelativistic reconnection studies. 



-17- 




Fig. 6. — Temporal evolution of the reconnection electric field (reconnection rate) at the 
X-point; (a) the raw reconnection rate r(t) = Ey/B^ and (h) the normalized reconnection 
rate f{t) = cEy/[cA,in'Bx,in']- The reference run L3 is presented in thick lines. The dotted 
lines contain negative mass density. 
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Figure [6] presents the reconnection rate r{t) and the normahzed reconnection rate f(t) 
for simulation runs L1-L9 in Table [1} Generally speaking, the lower-a runs L1-L5 last 
relatively long time. Their time duration t/xc ~ 200-300 is related to the plasmoid collisions. 
Therefore, the system has enough time to evolve to the Petschek-type reconnection, and we 
recognize Petschek-type current layers in these runs. We will visit the physical property of 
typical low-(T run (run LI) later in this section. Run L4 is the cold inflow counterpart of run 
L3; it uses the same parameters as run L3, except for the upstream plasma pressure. 

The higher-o" runs L6-L9 become unstable, and they stop before t/rc < 100. The numer- 
ical problem occurs around the plasmoids in the reconnection outflow front. As discussed 
in Section 3.1, there are discontinuities both in the upstream and the downstream of the 
magnetic pileup region. Since our numerical scheme (Lax-Wendroff scheme) is not ideal for 
shocks, we suffer from numerical noise at these discontinuities. Since the magnetic energy 
dominates the plasma energy in these runs, even small noises in the electromagnetic fields 
become crucial to fluid properties, and then the physically valid solution often collapses. 
In the dotted line region, we continue simulations even though we observe small negative 
mass in the edge of the plasma outflow, until our equation solver fails to find the mathe- 
matical solution. Since the numerical error occurs near the plasmoid, we think they show 
the right evolution for a while (20 ~ 30rc), until the unphysical information comes back to 
the X-point. We find that the normalized rate f{t) (Figure [6|)) is a better measure of the 
reconnection evolution, because it looks reasonably flat in higher-cx runs. However, we will 
only consider the times prior to the occurrence of negative density. 

As a general trend, we find that the reconnection rate becomes higher as the inflow 
density goes down, or as the parameter as^m increases. Figure [7] also shows the maximum 
reconnection rate f(t) in the simulation runs, as a function of the initial upstream a^^m 
parameter. In the limit of a^^in < 1, the reconnection rate is asymptotic to ~0.1. This 
is consistent with many studies on the nonrelativistic Petschek reconnection, whose the 
reconnection rate is known to be ~0.1. On the other hand, the rate constantly increases as 
the parameter a^^in increases. It is striking that the reconnection rate becomes closer to ~1, 
because the rate of one is the upper limit of magnetic dissipation. 

We briefly visit the global properties of low-cr runs. Top two panels in Figure |8] presents 
the late time snapshots at t/xc = 295 in run LI. Compared with the other higher-o" runs, 
the system evolution is rather slower due to the slow reconnection outflow. The typical 
outflow speed ~0.5c (0.57c at maximum) is consistent with the original upstream Alfven 
speed of 0.594c. In the current profile (Figure [8|)), we find Petschek-type current layers and 
the angle between current layers look wider than the reference run L3. Another current layer 
surrounding the plasmoid is very clear, too. These signatures are well observed in plasmoid 
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Fig. 7. — Dependence of the maximum reconnection rate r{t), as a function of the initial 
upstream parameter a^^in- 
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Fig. 8. — Top two panels: large-scale snapshots of run LI at t/rc = 295. (a) The x- 
component of the plasma 4-velocity Ux/c, and {h) the out-of-plane current jy/jo- Bottom 
two panels: snapshots of run L8 at t/Tc = 90. (c) The x-component of the inflow plasma 
4-velocity u^/c, and (d) the out-of-plane current jy/jo- 
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in nonrelativistic ion-electron plasmas. Unfortunately, we do not obtain long-term steady 
profile after the boundary collision, because the run stops immediately after this stage. 

The bottom two panels in Figure [8] show snapshots of the second most extreme case, 
run L8. In the outflow profile (Figure |8]c), we see that the outflow channel is narrower than 
the slower counterparts. At the edge of the Sweet-Parker outflow jets, the outflow 4- velocity 
becomes further relativistic, Ux/c ~ ±8.6, and the maximum Lorentz factor in the system 
is up to ~9. The current structure remains in a single current (Figure |8]d) at least at this 
stage. In the very thin current layer, there are small seeds of secondary tearing islands (e.g. 
a bright spot at x/L ~ —20 in the current profile; Figure |8](i). 



(a)Ey t= 80.0 xio° 
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Fig. 9. — Snapshots of run L9 at I/tc = 80. (a) The reconnection electric field Ey/ Bq. The 
white contour line indicates the region where the Lorentz invariant {E"^ — B"^) is positive, (b) 
The z-component of inflow 4- velocity Uz/c. 
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Furthermore, the most extreme case (run L9) shows an interesting evolution. Panels 
in Figure [9] show characteristic properties at t/r^ = 80, just before we meet an unphysi- 
cal solution at t/xc = 80.8. Figure |9]a shows the reconnection electric field. Importantly, 
its amplitude is Ey/Bo ~ 0.35-0.5. We indicate the "electric-dominant" region where the 
Lorentz invariant {E'^ — B"^) is positive with the white line. We generally observe such 
an electric-dominant region at the closer vicinity of the X-point, because the reconnection 
electric field Ey remains finite, while the magnetic field \B\ becomes zero at the X-point. 
However, such an electric-dominant region is usually confined in a very narrow region of the 
center of the reconnecting current sheet. For example, in run L8, such a region is very thin 
around the neutral plane, —0.5 < z/L < 0.5. However, in run L9, the electric field Ey be- 
comes so strong that it even dominates the magnetic field in a relatively large spatial region 
of —24 < x/ L < 24, —6 < z/ L < Q. In response to a strong electric field, we also find a 
super fast reconnection infiow (Figure |9]6). The maximum momentum is up to \uz\/c ^ 6.7, 
and the maximum infiow velocity is up to |f^|/c ~ 0.936. Also, as seen in Figure |9f , the re- 
connecting current layer becomes thicker 1-2L, while in other cases the central current layer 
always becomes thin < 0.5L. Since the plasma temperature becomes hot 5-lOmc^ along the 
current sheet, the kinetic scale increases by a factor of 5-10 and then it is comparable to the 
initial sheet thickness L. Therefore, we may have to consider kinetic effects beyond the fluid 
approximation. Indeed, an effective resistivity based on the kinetic effects is a long-standing 
problem in reconnection physics (e.g.. 

Regarding the energy conversion rates, we noticed that the magnetic pileup regions 
are also important in the relativistic runs. As a increases, the pileup flelds become more 
strong, and then more energy is delivered to the downstream Harris sheet plasmas there. 
On the other hand, in the current sheets and in the Petschek-type current layers, the energy 
conversion rate seems to be proportional to the reconnection rate. However, unfortunately, 
we do not have sufficient simulation results to discuss energy conversion in the high-a regime, 
which is of strong astrophysical interest. 



Hesse et al. (Tffl) 



3.3. Petschek-type current layer 



One of the most characteristic features of the late-time evolution of reconnection is 
the bifurcated Petschek-type current layers. We observe such current layers in runs L1-L5. 
In this section, we study how these current layers are influenced by the upstream energy 
composition CTe^in- In the relativistic Petschek reconnection, Lyubarsky (2005) examined the 
RMHD jump conditions across the slow shocks, and he found that the slow-shock angle 
becomes narrow when am in ^ 1- We examine our simulation results based on a similar 
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theory. Since Lyubarsky (2005)'s original work employs single-fluid RMHD model and it 
neglects the inflow plasma pressure, flrst, we construct complete jump conditions which 
contains both two-fluid effects and the inflow pressure. 




Fig. 10. — Rotated coordinate for the Petschek current layer. The thick black line stands 
for the current layer. The two angles 9^ and 9^ are relevant to the current layer and the 
upstream magnetic fleld line, respectively. The angle 9'^ is the field line angle in the rotated 
frame. 



Let us consider a rotated coordinate based on the Petschek current layer (or the slow 



shock surface in Lyubarsky (2005)). The new x' z' coordinate is tilted from the simulation 



coordinate xz by the angle of 9c as shown in Figure 10 The angles 9m and 9'^ are the 



upstream field line angles from the simulation frame and the rotated frame, respectively. 
Since the electric field is almost uniform over these regions in our simulation, we assume 
that the electric field Ey is constant. The relativistic stable conditions across the current 
layer are as follows: 



47r 
+ 2p + 
2wUx'Uzi Bx'Bz' 









(22) 
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Vx'Bz' — Vz'Bx' 

nuz' 



-- 
= 
0, 



(25) 
(26) 
(27) 



where the brackets stand for the jump condition in the z' direction. We employ the as- 
sumption of Vx'u = 0, Bx'd = 0, where u and d denote the upstream and the downstream 
properties. We confirmed that these assumption are fair, especially B^'d = 0. Then, equation 



25 yields 



Vz'u j3 Vx'd J-, 
^x'u — -Dz'- 



From equations 22 and 28 



+ -^J'^^'u = 2-fdWdUz'd- 
Bx'uBz' 2WdUx'dUz'd 



From equation 24 
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Eliminating Bz' with equation 28 
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From equations 29 and 31 , we obtain 
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^ c ' 2^lw^ + BlJAn 1 + cos2 
where we set a^^u = Bl/[An{2'ylwu)]- We also obtain 
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(28) 
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(30) 

(31) 
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(33) 



^ c ^ 2^lw^ + BlJAn l + a,,^cos^ej 
Then, we discuss the angles in the limit of as^u ^ 1- Approximating Vz'u = — ctan6'^ and 



Wd = Apd, equations 23 and 29 can be modified as follows: 
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We immediately obtain 
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It is reasonable that the outflow Lorentz factor is similar to that of upstream Alfven speed, 
lyr^FcT^. Considering that the outflow travels toward the +x direction, we find 

9, ~ 9'J{2a,,u). (39) 

This means that the Petschek outflow channel becomes narrower and narrower, as the up- 
stream flow is more and more magnetically dominated. 

In our simulation, we observe the Petschek-type current structures in runs L1-L5. In 
the other runs, as discussed, we could not solve the late-time evolution because of the 
numerical problems. In runs L1-L5, we measured the angle of the Petschek current layers in 
the following way. Near the current layer, we assume the inclined coordinate assuming an 



arbitrary angle 6(. like Figure 10 Then, across the current layer, we look at the relativistic 



jump conditions across the z' direction (eqs. 22 27). 



Figure 11 shows one example, physical properties across the current layers aX x/ L = 100 
at t/xc = 400, as indicated by the white line in Figure [2]d. In this case, the oblique frame 
properties are calculated by using an angle 9c = 0.125, and the opposite rotation is applied to 
the properties of the lower half and the upper half. We note that the neutral plane is slightly 
off-center (z/L ~ —0.15) in this very late stage because of the open boundary conditions. In 
the dense plasma region between the two current peaks, we observe fast reconnection outflow 



Vx' ~ 0.9c (Figure 11). We also observe noises in the properties near the center and the flux 
properties in the current layers; however, we think that they are sufficient for the purpose 
of this study. 

Varying 9c with /\9c = 0.025, we find out the best angle, which minimizes the variation 
of the above variables. Among them, the energy flux and the tangential momentum flux 



(Equations. 22 and 24 Figure 11c) in the outflow region and in the current layers are very 
sensitive, and so they give a reasonable estimate of 9c. We can also confirmed that the 
obtained angles are consistent with the topological structure, because the distance between 
the current peaks is ~2.4 and the location is x/L = 100. We repeat this procedure at various 
points along the well-developed current layers, where the structure is not influenced by the 
backward plasma flow around the plasmoids. Repeating the analyses at various time steps, 
we obtain the typical 9c angle for the specific run. 



Figure [12] compares the obtained angles by the above analysis in runs L1-L5. The 
dashed hne shows the current layer angle 9c- The typical field line angle 9m is also measured 
in the upstream side of the current layers, and they are presented in the solid line. We find 
that the angle 9c becomes narrower as the inflow parameter as^in increases. On the other 
hand, the field line angle 9^, shows the opposite trend. Considering that the reconnection 
rate increases as a^^in increases, it is quite reasonable that 9^ increases. We expect that 
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Fig. 11. — Physical properties across the Petschek-type current layers aX x/L = 100 at 
I/tc = 400. (a) Normalized plasma density 'jn/riQ, tangential plasma velocity v^'/c, normal 
plasma velocity v^'/c, (6) the out-of-plane electric current jy/jo, tangential magnetic field 

total pressure (eq. |23j normalized by 



Bx'/Bq, normal magnetic field B^i/Bq (eq. 27), (c) 
-Bq/Stt), and energy fiow (eq. 22 normalized by cBq/Sh) 
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Fig. 12. — Petschek angle analysis. The magnetic field line angle 9m {solid line), the current 
layer angle 9c {dashed line), and the estimated magnetic field line angle 9'^/2ae^in {dotted 
line) are presented as a function of the initial parameter (T£^j„. 
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the angle is eventually asymptotic to 9m ~ vr/4. The dotted line in Figure 12 shows the 
theoretical angle. It is estimated by substituting a^^u ~ o"£,jn in equation 39 In the present 
parameter range, we find an excellent agreement between two shock angles Oc and O'^l (2o"e,m) 
(the dotted line in Figure 12). Although we discuss stable current layers in mildly relativistic 
runs, we expect that the theory shows good agreement in the higher-o" regime, where the 
theory was originally designed. 



3.4. Uniform resistivity case 

In order to study the role of the resistivity, we also carried out another simulation run 
with a uniform resistivity (run U3 in Table [T]). The parameters are the same as those of 
run L3, but the resistivity is uniformly set. Its effective Reynolds number is Rm = 3000. 
Compared with run L3, the system evolves slower primary due to the low resistivity at the 



reconnecting X-point. The top three panels in Figure 13 present the late-time evolution of 



run U3, at I/tc = 200 and 300. Figure 13 i shows the properties along the outfiow line at 
t/Tc = 300. At t/Tc = 200, the reconnection outfiow is still only half way to the boundaries. 
The reconnecting current sheet contains several secondary structures. We think that this is 
due to the slower evolution of the system. There is sufiicient time for secondary structures 
to grow. The biggest plasmoids reach the boundaries around t/T^ = 300. Now we observe 
a formation of multiple big islands inside the reconnecting current sheet. As wee see in 
the profiles in Figure [T3]d, multiple magnetic reconnections take place and expel outfiows 



between these islands. The outfiow 4- velocity reaches u^/c ~ 3 at various local points, and 



the global fiow speed seems to be m^^/c ~ 1-2. The density spikes in Figure 13 i are identical 
to the 0-points, magnetic nulls at the center of plasmoids. Although the out-of-plane fiow is 
very small, Uy/c <^ 1, these high-density plasmas carry the electric current inside the 0-type 



regions (Figure 13 c). On the other hand, around several regions between the islands, we 



see that the out-of-plane 4- velocity is enhanced, Uy/c ~ 1 or 1.5 (Figure 13d). They are 
related to thin current sheets between plasmoid islands. The plasma temperature is typically 
p/nmc^ ~ 2 in the outfiow region, and it becomes very high p/nmc^ ~ 4-5 around the O- 
points. The simulation continues until I/tc ^ 345 shortly after the plasmoids completely 
went through the boundaries. 

When the plasmoid islands appear, its typical timescale seems to be tens of Tc, and it 
is faster than an estimated timescale of the resistive tearing mode, -R^^ or hE^l^ ~ (9(10^). 
We think that the island formation is enhanced by the two-fiuid effect, which was introduced 



in our simulation. Since our Ohm's law (eq. 18) contains the fluid inertial term dt{hpU. 



bp Uipy J , 



the tearing mode can grow more explosively than the classical resistive MHD case. If we use 
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Fig. 13. — Large-scale snapshots of run U3: (a) the x-component of the plasma 4-velocity 
Ux/c at t/Tc = 200, (6) the x-component of the plasma 4-velocity Ux/c at t/rc = 300, and (c) 
the out-of-plane current jy/jo at t/Tc = 300. The black lines show magnetic field lines, (d) 
The outflow 4-velocity Ux/c, the out-of-plane 4-velocity Uy/c, and the plasma density 'jn/rio 
at t/rc = 300, along the right half of the outflow line {z = 0). 
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the specific condition oi Uy ~ 1, the timescale of the relativistic collisionless tearing mode 
( [Zelenyi fc Krasnoselskikh|[l979| is C(10), too. 

Although the system evolution is slower than the reference run L3, we note that the 
reconnection still remains fast, at least during this simulation run, and it may be related 
to island formation. If we discuss the global structure by filtering out the local plasmoid 
islands, the average plasma inflow speed is Vy ~ 0.1c, and the reconnection electric field 
is Ey ~ 0.1i?o- This may be an interesting hint to discuss the problem of a fast magnetic 
reconnection. 



4. DISCUSSION 



First, let us briefly compare our results with the one-fluid work by Watanabe fc Yokoyama 



(2006) (here referred as WY06). They employ an relativistic Ohm's law 



E 



V 



X B = {ri/i)j, 



(40) 



where 7' is the Lorentz factor of the one-fluid MHD motion. Our Ohm's equation (eq. 18) 
differs in the following two ways. First, since our equation contains the fluid inertial term, by 
definition our model describes better physics. Second, we do not consider the factor of I/7'. 
However, we consider finite resistivity only near the X-point, where 7' is close to the unity. 
The fastest run in WY06 is directly equivalent to our reference run L3. The reconnection 
geometry looks similar. However, we find various minor differences. The maximum outflow 
4-velocity is u^/c ~ 2.2-2.3 in WY06, while we often observe faster value {u^/c > 3) (e.g.. 
Figure [3^). This is quite probably due to the two-fluid effect and the grid condition. Since 



we also deal with the out-of-plane motion Uy, the Lorentz factor can be larger even when it 
contains a contribution from Vy. Also, WY06 employed nonuniform grids, and then physical 
quantities in the distant outflow region are often averaged in the larger computational cells. 
Regarding the structure of the Petschek-type reconnection, WY06 implied that their shock 



angle becomes narrower in the relativistic regime like Lyubarsky (2005) predicted. Taking 



the inflow pressure into account, we clarified that a modified Lyubarsky (2005) theory well 



explain the simulation results (Figure 12). Our angle is narrower than that of WY06. This is 
probably because the isotropic plasma pressure is usually overemphasized in one-fluid model, 
and because the amplitude of the effective resistivity may be different. 

We obtained several implications for relativistic reconnection models. Although we do 
not obtain a steady state Sweet-Parker reconnection, in all our runs, the early evolutions 
of main reconnection runs will be good hints to understand relativistic Sweet-Parker recon- 
nection. For example, in run L3, along the outflow line, we found that plasma temperature 
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becomes very high (e.g., Figure |3]e). Consequently, the relativistic enthalpy substantially 
increases. Aw ~ 4nAT > nmc^ in the outflow region. Therefore, the relativistic recon- 



nection model cannot neglect the relativistic gas pressure in the outflow region (Lyubarsky 



2005 Zenitani & Hesse 2008b). In the Petschek reconnection regime, we observed that the 



bifurcated current layers and their angle becomes narrower and narrower as the inflow be- 
comes more and more magnetically dominated. With minor modification, this trend is well 



consistent with the argument proposed by Lyubarsky (2005). 



The ultrarelativistic limit of as^in 3> 1 is of strong astrophysical interests. Our parameter 
study suggested that the reconnection rate is asymptotic to its upper hmit of ~1 in that 
regime. Indeed, in the most extreme case (run L9), we observe fast reconnection with super 
fast inflow. We expect that reconnection is super fast in the high-a regime. Such a fast 
reconnection rate implies that the separatrix angle will be wide open — asymptotic to 45° in 
the steady stage. In fact, our mildly relativistic runs show the magnetic field line angle 9m 



constantly increases (e.g.. Figure 12) in the Petschek-type steady regime. In a sense this 
is reasonable, because there are less current carrier in such a regime. When the separatrix 
becomes open, the field reversal current for the reconnected fields ±5^ partially cancel the 
field reversal current for the antiparallel fields iiB^, therefore the system needs less electric 
current. We do not know whether or not the bifurcated Petschek-type solution exists in the 
high-cT regime. Since the current layer becomes too flat and the central diffusion region tends 
to expand, the reconnection current sheet may remain in a single thick current layer for a 
long time. 

We think that an important feature of the reconnection in the high-o" regime is the 
shortage of the current carrier. As the authors discussed through PIC simulation and the 



two- fluid theory (Zenitani & Hesse 2008b), when reconnection environment is magnetically 



dominated and runs out of current carriers, the displacement current induces the strong 
electric field. It leads to a faster reconnection rate and the expansion of the central Sweet- 
Parker region. From MHD viewpoint, it means the enhancement of the effective resistivity; 
however, we note that the conventional single-fluid RMHD simulations have no explicit upper 
limit of plasma currents. Due to the enhancement of the reconnection field, we find a large 
electric-dominated region in run L9, where the field is electrically dominated, {E'^ — B^) > 0, 
around the X-type region. Plasmas are no longer magnetized there, and then powerful 
DC acceleration will occur ( Zenitani fc Hoshino||2001 ). Therefore, we expect that the high-cx 
reconnection is a favorable source of nonthermal particles acceleration. Long-term evolution, 
theoretical modeling, and particle acceleration in the high-cx regime will be left for future 
work. 



In this work, we mainly use the energy-based magnetization parameter ae,in, because it 
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seems to be a better measure of reconnection property than the conventional magnetization 
parameter cx^.m- In fact, as long as we surveyed (Sections 3.2 and 3.3), run L4 with lower 
pressure resembles run L5 rather than run L3. However, we still observe minor differences, 
and the running out of the current carrier should be controlled by am- So, we conclude that 
both two magnetization parameters as,in and cr^.m are important. 

We also demonstrated that the spatial profile of the resistivity has great influence on the 
system evolution. As recognized in many works in the nonrelativistic regime, the spatially 
localized resistivity leads to the Petschek-like reconnection with bifurcated current layers 



(e.g., Ugai & Tsuda (1977); Scholer (1989)). On the other hand, the uniform resistivity 
case exhibits a single current layer with secondary islands. Considering that the governing 
equations are almost the same outside the localized resistivity point, it is impressive to 
see such a contrast. The dependence on the amplitude of the resistivity r^e//, the spatial 
profile, and the other physical models, will be left for future work. Regarding the outflow 
structure, at present PIC simulations of the relativistic magnetic reconnection exhibit a 
laminar outflow without islands or with small minor islands in a main reconnection region 



(Zenitani & Hoshino 


2001, 


2007 


Jaroschek et al. 


2004 


Zenitani & Hesse 


2008a|b 


). On the 


other hand, in the nonrelativistic regime. 


Daughton & Karimabadi 


(2007 


) demonstrated a 



very interesting result by using a large-scale PIC simulation. They showed that reconnection 
outflow is highly influenced by a continuous formation of secondary islands, and its global 



picture looks similar to our uniform resistivity run (Figure 13b). We do not know whether 
the relativistic reconnection is influenced by such continuous island formation, it is worth 
investigating by using a larger relativistic PIC simulation. We also observe small islands in 
higher-o" runs, and so the island formation may also be controlled by the upstream parameters 
(c"m,i?i and 0"^ j„). 

On the viewpoint of numerical accuracy, we confirmed that our primitive variable solver 



is sufficiently accurate. In Appendix A (Figure 14), the numerical error of our solver is 
presented. In our range of interest {\u\/c ^ 10^ and p/nmc^ < 10), the relative errors in the 
restored primitive variables are very small, ~ O{10~^^). Therefore, the worst estimate of the 
accumulated error would be still negligible, 10^ x 10"^^ ■ (400rc/At) < 10"^ On the other 
hand, in order to further explore the higher-o" conditions, we have to improve the numerical 
scheme. At present, the modified Lax-Wendroff scheme seems to be the bottle neck. It is 
not ideal to describe shocks, while the discontinuities around the magnetic pileup regions 
are always difficult to solve in a nonsteady stage of reconnection. We plan to employ a more 



stable scheme such as HLL schemes (Mizuno et al. 2006 Mignone et al. 2009) in order to 



study long-term evolution in the high-a regime. 



Finally, let us discuss potential targets beyond this work. A straightforward extension 
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will be magnetic reconnection with the out-of-plane magnetic field [By or the "guide field"). 
In is already known that PIC simulations with the guide field exhibit global charge separation 
in a reconnection region, and so the neutral one-fiuid approximation already breaks down 
(Zenitani & Hoshino 2008). Therefore we have to solve the positron and electron evolution 
separately without the symmetric assumption. In the Petschek regime, Lyubarsky (2005) 
claimed that (1) the relativistic magnetic reconnection involves rotational discontinuities as 
well as slow shocks, and that (2) the compressed guide field fiux inside the outfiow channel is 
the main energy carrier in the magnetically dominated regime. These properties are worth 
checking in future simulations. In three dimensions, it is known that the reconnection current 
sheet is unstable to the relativistic drift kink instability, which arises from the counter- 
streaming two-fiuid motion of positron fiuids and electron fiuids ( Zenitani fc Hoshino||2005a 



2007 Pritchett et al. 1996 Daughton 1999). These three-dimensional evolutions should be 



carefully compared with PIC simulations, so that we can study the larger problems such as 
magnetar fiares and global pulsar magnetospheres by using a relativistic two-fiuid model. Of 
course, it is critically important to establish an improved (theoretical or empirical) resistivity 
model, which highly affects the system evolution. 



5. SUMMARY 

We carried out relativistic two-fluid MHD simulation of a magnetic reconnection in an 
electron-positron pair plasma. The interspecies friction term works as an effective resistiv- 
ity, and then we successfully demonstrated the large-scale evolution of relativistic magnetic 
reconnection. The system evolves from a Sweet-Parker-like fast reconnection to a Petschek- 
like reconnection with bifurcated current layers. Open boundary conditions enable us to 
observe long-term evolutions, and we flnd a Petschek structure, which is quite stable. As 



Lyubarsky (2005) predicted, the current layer angle becomes substantially narrower when 
the reconnection inflow is more magnetically dominated. Meanwhile, we flnd that the recon- 
nection rate goes up to ~1 in extreme cases, which implies that efficient particle acceleration 
occurs in the electric-dominated region. In addition, we demonstrate that the system evolu- 
tion is controlled by the resistivity model. We emphasize that the large-scale reconnection 
problems are investigated with a two-fluid RMHD model. Beyond the single-fluid RMHD 
approximation, multifluid models will be good alternatives to study astrophysical plasma 
problems which involve magnetic dissipation. 
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Table 1: List of Simulation runs 



Name 


Domain Size 


Grid Points 


CAt/Ag 




Ptn/Po 








S3 


80 


X 


40 


1200 


X 


600 


0.3 


0.1 


1.0 


20 


A 


0.894 


M3 


240 


X 


120 


3600 


X 


1800 


0.3 


0.1 


1.0 


20 


A 


0.894 


LI 


240 


X 


120 


4^ 


^00 


X 


2400 


0.2 


1.0 


1.0 


2 


0.4 


0.535 


L2 


240 


X 


120 


4^ 


^00 


X 


2400 


0.2 


0.3 


1.0 


6.67 


1.33 


0.816 


L3 


240 


X 


120 


4^ 


^00 


X 


2400 


0.2 


0.1 


1.0 


20 


4 


0.894 


L4 


240 


X 


120 


4^ 


^00 


X 


2400 


0.2 


0.1 


0.3 


20 


9.1 


0.953 


L5 


240 


X 


120 


Ai 


^00 


X 


2400 


0.2 


0.05 


1.0 


40 


8 


0.949 


L6 


240 


X 


120 




^00 


X 


2400 


0.2 


0.03 


1.0 


66.7 


13.3 


0.964 


L7 


240 


X 


120 


Ai 


^00 


X 


2400 


0.2 


0.02 


1.0 


100 


20 


0.976 


L8 


240 


X 


120 


Ai 


^00 


X 


2400 


0.2 


0.01 


1.0 


200 


40 


0.988 


L9 


240 


X 


120 


4^ 


^00 


X 


2400 


0.2 


0.005 


1.0 


400 


80 


0.994 


U3 


240 


X 


120 


4^ 


^00 


X 


2400 


0.2 


0.1 


1.0 


20 


4 


0.894 


XL3 


240 


X 


120 


9600 


X 


4800 


0.2 


0.1 


1.0 


20 


4 


0.894 



*uniform resistiv 



A. Obtaining primitive variables 



The rest energy density the momentum density m, and the energy density £ are 
related to the primitive variables in the following way: 



D = ^nmc 
m = [j^e + p)u]/c^ 
£ = 7^(e + ]9) - p. 



(Al) 
(A2) 
(A3) 



For convenience, we introduce M = \m\c and u = \u\/c. We consider the case of M > 0, 
because we immediately know u = when M = 0. Approximating the enthalpy {e + p) = 
nmc^ + Gp by G = r/(r — 1), equations A2 and A3 become 



M 
£ 



'j{nmc^ + Gp)u = (D + 'jGp)u 
•y'^^nmc^ + Gp) — p = 7-D + (7^G — l)p 



From equations A4 and A5 , we can eliminate p in the following way: 

{-f^G - 1)M - -fGu£ = {-f^G - l)Du - -f^GuD 



(A4) 
(A5) 



(A6) 



- 35 - 



-iGu£ = (72^ - 1)M + Du. 
Squaring equation A7 and substituting 7^ = 1 + -u^, we obtain 



(A7) 



G'{S' - M')u^ - 2GMDu^ + - 2GM\G - 1 



2D{G- l)M u- {G-l)M 



We solve this equation to obtain the physically valid solution. Other primitive variables are 
easily obtained by using the solution u. When the first coefficient {E"^ — M^) is negative 
we immediately stop the simulation, because such situation is physically invalid. We also 
checked the other conditions D > Q and {S — D) > 0. 

The behavior of our primitive variable solver is characterized by two parameters, the 
relativistic bulk flow u and the relativistic temperature p/nmc^. We benchmarked the nu- 



merical accuracy of our solver, and Figure 14 shows the results as a function of the two 



parameters. We can see that the quartic solution u is accurate even in the ultrarelativistic 
regime of m ~ 10^. The pressure p is least reliable in the limit of p <^ nm(? and 1. This 
is because the pressure is enclosed in the enthalpy term w = e + p = n[mc^ + G{p/nmc^)\, 
but the fluid macro properties are insensitive to p in such cases. The error in n shows the 
same trend as that of u. 



(a) Error : u 

4 1 



■5 
1^ 



10-1° 10 = 



-2 2 4 

iogio(u/c) 



10° 

10-= 
10-10 

10-' = 



(b) Error : p 

44 



c 

Q. 



10-1= 10-' 



.J 



-4-2 2 

logio(u/c) 



Fig. 14. — Numerical accuracy of our primitive variable solver. Relative errors in the 
restored (a) u and {h) p are presented as a function of u and p/nmc^. 



B. Brown Method 



We solve the quartic equations by using a simplified version of the Brown method 



(|Nunohiro & Hirano||2003|). Consider the following quartic equation: 





+ cau^ + C2U^ + ciu + Co 



(Bl) 
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where c„ is the nth-order coefficient. When we apply eq. |A8| to eq. |B1[ we find that the 



coefficients c„ are all real, and that Cq < 0. In addition to two solutions (ui,M2), eq. Bl has 



one real negative solution {u^ < 0) and one real positive solution (0 < U4). The first two 
(mi,M2) are usually complex conjugates, and the real positive one U4 is the physically valid 
solution that we are looking for. Let us consider the following two equations: 



{u - Ui){u - U2) 
{u - u^){u - U4) 



u + aiu + 
+ hiu + 60. 



(B2) 
(B3) 



Examining the properties of hn and c„, we know that a„ G M and that ao > and 60 < 0. 



Therefore the quartic equation (eq. A8) can be decomposed into two quadratic equations 

and follows: 



with real coefficients. The relations between a^, 6„, 



C3 = ai + hi 

C2 = aibi + ao + bo 

Cl = a^bi + aibo 

Co = o-o^o- 



(B4) 



We further define v :- 
equations 



(fli — bi),w := (ao + bo), X := (ao — bo). Using v,w,x, we rewrite these 




c| — f ^ + Aw 
2C3W — 2xv 
w"^ — x'^. 



We obtain the following relation: 



2 2 

X V 



(csu; - 2cif = {w'^ - Aco){Aw - 4c2 + C3) 



(B5) 



(B6) 



and then 



w' 



- C2W^ + (C1C3 - Aco)w + [co(4c2 - C3) - cj 



(B7) 



We solve this third-order equation to obtain w, paying attention to the numerical accuracy 



(Nunohiro et al. 1996). Usually, we obtain two complex solutions and the one real solution 



in Equation B7, and so we employ the sole real solution w. When we find three real solutions 



in Equation B7, we may have multiple choices for u, because Equation B2 also has two real 
solutions instead of complex conjugates. 



By using this w, we obtain x = +\^w'^ — Aco- Here we choose a positive square root 
because we know that ao > 0, 60 < 0. The last variable v can be obtained from eq. B^ 
accordingly. By using w, x, v, and C3, we obtain a„ and 6„, and then we obtain M4 as a 



positive root of eq. B3 
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If we employed the standard Ferrari's method to solve the quartic equations, we should 
use Cardano's transformation, u := u' — C3/4, in order to eliminate the third-order coefficient 
C3. Then, the inverse transformation of u often causes the cancellation of significant digits 
when the solution is very small, m ~ 0. On the other hand, by using Brown method, we 
are not so influenced by the cancellation of signiflcant digits when the solution is very small, 
■u ~ 0. 
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